#### Setup R

# Clear terminal
cat("\014")

# Clear space
rm(list = ls())

# Load library packages
library(ggstatsplot)
library(RColorBrewer)
library(rio)
library(tidyverse)
library(broom)
library(ggplot2)
library(gridExtra)
library(hrbrthemes)
library(stringi)
library(tufte)
library(summarytools)
library(jtools)
library(ggthemes)
library(estimatr)
library(tidyverse)
library(purrr)
library(ggplot2)
library(gridExtra)
library(knitr)
library(modelsummary)
library(sjPlot)

#### Setup data 
dat <- read.csv("election-monitoring.csv") %>%
  mutate(time = as.numeric(Duration..in.seconds.))
dim(dat)
table(dat$T)

head(dat$StartDate)
tail(dat$EndDate)

cutoff <- median(dat$time, na.rm = T)/3

dat <- dat %>%
  filter(time >= cutoff) # dropping speeders per Lucid's guidance and metric
dim(dat)

# order dependent variable levels
dat$DV1 <- factor(dat$Q536, 
                  levels = c("Strongly disagree",
                             "Somewhat disagree",
                             "Neither agree nor disagree",
                             "Somewhat agree",
                             "Strongly agree"))

# create numerical variables
dat$DV1num <- as.numeric(dat$DV1)
dat$Q538 <- as.numeric(dat$Q538)
dat$probability <- dat$Q538
dat$Q538 <- NULL

table(dat$DV1num, dat$DV1)

prop.table(table(dat$DV1))

#construct indicator variables for treatment groups
dat$treatment <- dat$T
dat$treatment[dat$T == "C"] <- "Control"
dat$treatment[dat$T == "T1"] <- "Trump"
dat$treatment[dat$T == "T2"] <- "Biden"
dat$treatment[dat$T == "T3"] <- "IO"
levels(dat$treatment) <- c("Control", 
                           "Biden", 
                           "IO", 
                           "Trump")


dat$Trump <- as.numeric(dat$T == "T1")
dat$Biden <- as.numeric(dat$T == "T2")
dat$IO <- as.numeric(dat$T == "T3")

# recode party identification
dat$pid <- dat$Q16
dat$pid[dat$Q16 == "Independent"] <- NA
dat$pid[dat$Q16 == "Other"] <- NA
dat$pid <- as.factor(dat$Q16)

# setup colors
myColors <- rev(c("red4", "seagreen4", "black", "royalblue4"))
reds <- rev(c("red4", "red3", "red2", "red1"))
blues <- rev(c("royalblue4", "royalblue3", "royalblue2", "royalblue1"))

agreement <- c("Strongly agree",
               "Somewhat agree",
               "Neither agree nor disagree",
               "Somewhat disagree",
               "Strongly disagree")

#### Descriptives
dat%>% 
  drop_na(pid) %>%
  ggplot(aes(x = pid))+
  geom_bar(aes(y = (..count..)))+xlab("Party Identification")+ylab("Count")+geom_text(stat='count', aes(label=..count..), vjust=-1)+ylim(0, 1500) +
  theme_tufte(base_family = "GillSans")
ggsave("~/Dropbox/curtice-crabtree/plots/pid.png",
       height = 5,
       width = 5)

## Concerns about election fraud, intimidation, and violence
table(dat$Q33_1)
dat$violence<-factor(dat$Q33_1, levels=c("Don’t know", "Not at all concerned","Not very concerned","Somewhat concerned", "Very concerned"))
table(dat$violence)

table(dat$Q33_2)
dat$voter_suppression_intim<-factor(dat$Q33_2, levels=c("Don’t know", "Not at all concerned","Not very concerned","Somewhat concerned", "Very concerned"))

table(dat$Q33_3)
dat$election_rigging<-factor(dat$Q33_3, levels=c("Don’t know", "Not at all concerned","Not very concerned","Somewhat concerned", "Very concerned"))
table(dat$violence,dat$pid)
table(dat$Q33_4)
dat$not_free_fair<-factor(dat$Q33_4, levels=c("Don’t know", "Not at all concerned","Not very concerned","Somewhat concerned", "Very concerned"))

prop.table(table(dat$violence))
prop.table(table(dat$voter_suppression_intim))
prop.table(table(dat$election_rigging))
prop.table(table(dat$not_free_fair))

# crosstab concerns by partisanship
table(dat$violence,dat$pid)
table(dat$election_rigging,dat$pid)
table(dat$voter_suppression_intim,dat$pid)
table(dat$not_free_fair,dat$pid)

## Descriptive plots
ggplot(dat, aes(x = violence)) +  
  geom_bar(aes(y = (..count..)/sum(..count..)))

ggplot(dat, aes(x = election_rigging)) +  
  geom_bar(aes(y = (..count..)/sum(..count..)))

ggplot(dat, aes(x = voter_suppression_intim)) +  
  geom_bar(aes(y = (..count..)/sum(..count..)))

ggplot(dat, aes(x = not_free_fair)) +  
  geom_bar(aes(y = (..count..)/sum(..count..)))

## Concerns about election fraud, intimidation, and violence
table(dat$Q33_1)
dat$"Political Violence"<-factor(dat$Q33_1, levels=c("Don’t know", "Not at all concerned","Not very concerned","Somewhat concerned", "Very concerned"))
table(dat$violence)

table(dat$Q33_2)
dat$"Voter suppression/intimidation"<-factor(dat$Q33_2, levels=c("Don’t know", "Not at all concerned","Not very concerned","Somewhat concerned", "Very concerned"))

table(dat$Q33_3)
dat$"Election Rigging"<-factor(dat$Q33_3, levels=c("Don’t know", "Not at all concerned","Not very concerned","Somewhat concerned", "Very concerned"))
table(dat$violence,dat$pid)
table(dat$Q33_4)
dat$"Election will not be free or fair"<-factor(dat$Q33_4, levels=c("Don’t know", "Not at all concerned","Not very concerned","Somewhat concerned", "Very concerned"))


myvars <- c("Political Violence","Voter suppression/intimidation","Election Rigging","Election will not be free or fair")
dat1 <- dat[myvars]
dat2<-dat1 %>%  
  gather(key = items, value = answer) %>% 
  mutate(answer = factor(answer)) 
dat2<-data.frame(dat2)
dat2$answer<-factor(dat2$answer,labels=c("Don’t know", "Not at all concerned","Not very concerned","Somewhat concerned", "Very concerned"))

p1 <- dat2 %>% 
  drop_na(answer) %>%
  ggplot(aes(x = (items))) +
  geom_bar(aes(fill = answer), position = "fill")

p1 <- p1 + 
  coord_flip() + 
  scale_fill_brewer(palette = "Greys") + 
  labs(title = "Thinking ahead to the presidential election, how concerned are you about the following things?",
       subtitle = "All Respondents") +
  xlab("") + 
  ylab("Proportion") +  
  theme_tufte(base_family = "GillSans")
p1
ggsave("~/Dropbox/curtice-crabtree/plots/desc.png",
       p1,
       height = 5,
       width = 10)

#### First outcome measure - Agreement
treatments.vis <- ggstatsplot::ggbetweenstats(
  data = dat,
  x = treatment,
  y = DV1num,
  xlab = "Treatments",
  ylab = "Agreement",
  plot.type = "violin",
  type = "np",
  pairwise.display = "significant",
  centrality.type = "parametric",
  conf.level = 0.95,
  title = "",
  package = "ggsci",
  palette = "nrc_npg"
) + 
  geom_jitter(alpha = 0.05) +
  theme_tufte(base_family = "GillSans") +
  theme(legend.position = "none") +
  scale_colour_manual(values = myColors) +
  labs(subtitle = "Agreement: average treatment effects")
treatments.vis
ggsave("~/Dropbox/curtice-crabtree/plots/ate.png",
       height = 7,
       width = 9)

###### TEH by PID
summary(lm_robust(DV1num ~ treatment*pid, data = dat))

treatments.rep.vis <- ggstatsplot::ggbetweenstats(
  data = dat[dat$pid == "Republican", ],
  x = treatment,
  y = DV1num,
  xlab = "Treatments",
  ylab = "Agreement",
  plot.type = "violin",
  type = "np",
  pairwise.display = "significant",
  centrality.type = "parametric",
  conf.level = 0.95,
  title = "",
  package = "ggsci",
  palette = "nrc_npg"
) + 
  geom_jitter(alpha = 0.05) +
  theme_tufte(base_family = "GillSans") +
  theme(legend.position = "none") +
  scale_colour_manual(values = myColors) +
  labs(subtitle = "Agreement: Conditional average treatment effects (Republicans)")
treatments.rep.vis
ggsave("~/Dropbox/curtice-crabtree/plots/cate-rep.png",
       height = 7,
       width = 9)

treatments.dem.vis <- ggstatsplot::ggbetweenstats(
  data = dat[dat$pid == "Democrat", ],
  x = treatment,
  y = DV1num,
  xlab = "Treatments",
  ylab = "Agreement",
  plot.type = "violin",
  type = "np",
  pairwise.display = "significant",
  centrality.type = "parametric",
  conf.level = 0.95,
  title = "",
  package = "ggsci",
  palette = "nrc_npg"
) + 
  geom_jitter(alpha = 0.05) +
  theme_tufte(base_family = "GillSans") +
  theme(legend.position = "none") +
  scale_colour_manual(values = myColors) +
  labs(subtitle = "Agreement: Conditional average treatment effects (Democrats)")
treatments.dem.vis
ggsave("~/Dropbox/curtice-crabtree/plots/cate-dem.png",
       height = 7,
       width = 9)

treatments.ind.vis <- ggstatsplot::ggbetweenstats(
  data = dat[dat$pid == "Independent" |
               dat$pid == "Other", ],
  x = treatment,
  y = DV1num,
  xlab = "Treatments",
  ylab = "Agreement",
  plot.type = "violin",
  type = "np",
  pairwise.display = "significant",
  centrality.type = "parametric",
  conf.level = 0.95,
  title = "",
  package = "ggsci",
  palette = "nrc_npg"
) + 
  geom_jitter(alpha = 0.05) +
  theme_tufte(base_family = "GillSans") +
  theme(legend.position = "none") +
  scale_colour_manual(values = myColors) +
  labs(subtitle = "Agreement: Conditional average treatment effects (Independents/Others)")
treatments.ind.vis
ggsave("~/Dropbox/curtice-crabtree/plots/cate-ind.png",
       height = 7,
       width = 9)

g <- grid.arrange(treatments.dem.vis,
                  treatments.ind.vis,
                  treatments.rep.vis,
                  nrow = 3)
ggsave("~/Dropbox/curtice-crabtree/plots/cate-all.png",
       g,
       height = 14,
       width = 9)

#### Regression models for agreement outcome 
lmout <- lm_robust(DV1num ~ Trump + Biden + IO + Q21 + Q25 + Q23 + as.numeric(Q24), data = dat)
summary(lmout)

lmout.dem <- lm_robust(DV1num ~ Trump + Biden + IO + Q21 + Q25 + Q23 + as.numeric(Q24), data = dat[dat$pid == "Democrat", ])
summary(lmout.dem)

lmout.io <- lm_robust(DV1num ~ Trump + Biden + IO + Q21 + Q25 + Q23 + as.numeric(Q24), data = dat[dat$pid == "Independent" |
                                                                                                           dat$pid == "Other", ])
summary(lmout.io)

lmout.rep <- lm_robust(DV1num ~ Trump + Biden + IO + Q21 + Q25 + Q23 + as.numeric(Q24), data = dat[dat$pid == "Republican", ])
summary(lmout.rep)

terms <- names(coef(lmout))

plot_models(lmout, lmout.io, lmout.dem, lmout.rep,
                   rm.terms = terms[-c(2:4)],
                   show.values = T,
                   show.p = T,
                   legend.title = "Models",
                   m.labels = c("All", "Independents/Others", "Democrats", "Republicans"),
                   vline.color = "lightgrey",
                   p.adjust = "holm",
                   colors = "quadro",
                   digits = 3) +
  labs(title = "Agreement: Treatment effects") + 
  ylab("") + 
  xlab("Treatment Effects") + 
  theme_tufte(base_family = "GillSans") 
ggsave("~/Dropbox/curtice-crabtree/plots/ols.png",
       height = 7,
       width = 7)

cm <- c('Trump'    = 'Trump',
        'Biden'    = 'Biden',
        'IO' = 'IO')

modelsummary(list(lmout, lmout.io, lmout.dem, lmout.rep),
             output = "latex",
             coef_map = cm,
             stars = T)


dat$Party<-factor(dat$Q16,levels = c("Democrat","Republican","Independent","Other"))
party_subsets <- 
  dat %>%
  split(.$Party) %>%
  map( ~ lm_robust(DV1num ~ Trump + Biden + IO + Q21 + Q25 + Q23 + as.numeric(Q24), data = .)) %>%
  map(tidy) %>%
  bind_rows(.id = "Party")

#### Second outcome measure - margin
treatments.vis.mar <- ggstatsplot::ggbetweenstats(
  data = dat,
  x = treatment,
  y = probability,
  xlab = "Treatments",
  ylab = "Margin of victory",
  plot.type = "violin",
  type = "np",
  pairwise.display = "significant",
  centrality.type = "parametric",
  conf.level = 0.95,
  title = "",
  package = "ggsci",
  palette = "nrc_npg"
) + 
  geom_jitter(alpha = 0.05) +
  theme_tufte(base_family = "GillSans") +
  theme(legend.position = "none") +
  scale_colour_manual(values = myColors) +
  labs(subtitle = "Margin of victory: average treatment effects")
treatments.vis.mar
ggsave("~/Dropbox/curtice-crabtree/plots/ate-mar.png",
       height = 7,
       width = 9)

treatments.rep.vis.mar <- ggstatsplot::ggbetweenstats(
  data = dat[dat$pid == "Republican", ],
  x = treatment,
  y = probability,
  xlab = "Treatments",
  ylab = "Margin of victory",
  plot.type = "violin",
  type = "np",
  pairwise.display = "significant",
  centrality.type = "parametric",
  conf.level = 0.95,
  title = "",
  package = "ggsci",
  palette = "nrc_npg"
) + 
  geom_jitter(alpha = 0.05) +
  theme_tufte(base_family = "GillSans") +
  theme(legend.position = "none") +
  scale_colour_manual(values = myColors) +
  labs(subtitle = "Margin of victory: Conditional average treatment effects (Republicans)")
treatments.rep.vis.mar
ggsave("~/Dropbox/curtice-crabtree/plots/cate-rep-mar.png",
       height = 7,
       width = 9)

treatments.dem.vis.mar <- ggstatsplot::ggbetweenstats(
  data = dat[dat$pid == "Democrat", ],
  x = treatment,
  y = probability,
  xlab = "Treatments",
  ylab = "Margin of victory",
  plot.type = "violin",
  type = "np",
  pairwise.display = "significant",
  centrality.type = "parametric",
  conf.level = 0.95,
  title = "",
  package = "ggsci",
  palette = "nrc_npg"
) + 
  geom_jitter(alpha = 0.05) +
  theme_tufte(base_family = "GillSans") +
  theme(legend.position = "none") +
  scale_colour_manual(values = myColors) +
  labs(subtitle = "Margin of victory: Conditional average treatment effects (Democrats)")
treatments.dem.vis.mar
ggsave("~/Dropbox/curtice-crabtree/plots/cate-dem-mar.png",
       height = 7,
       width = 9)

treatments.ind.vis.mar <- ggstatsplot::ggbetweenstats(
  data = dat[dat$pid == "Independent" |
               dat$pid == "Other", ],
  x = treatment,
  y = probability,
  xlab = "Treatments",
  ylab = "Margin of victory",
  plot.type = "violin",
  type = "np",
  pairwise.display = "significant",
  centrality.type = "parametric",
  conf.level = 0.95,
  title = "",
  package = "ggsci",
  palette = "nrc_npg"
) + 
  geom_jitter(alpha = 0.05) +
  theme_tufte(base_family = "GillSans") +
  theme(legend.position = "none") +
  scale_colour_manual(values = myColors) +
  labs(subtitle = "Margin of victory: Conditional average treatment effects (Independents and Others)")
treatments.ind.vis.mar
ggsave("~/Dropbox/curtice-crabtree/plots/cate-dem-mar.png",
       height = 7,
       width = 9)

g <- grid.arrange(treatments.dem.vis.mar,
                  treatments.ind.vis.mar,
                  treatments.rep.vis.mar,
                  nrow = 3)
ggsave("~/Dropbox/curtice-crabtree/plots/cate-all-pr.png",
       g,
       height = 14,
       width = 9)

### Regression models for margin of victory outcome
lmout.probability <- lm_robust(probability ~ Trump + Biden + IO + Q21 + Q25 + Q23 + as.numeric(Q24), data = dat)
summary(lmout.probability)

lmout.dem.probability <- lm_robust(probability ~ Trump + Biden + IO + Q21 + Q25 + Q23 + as.numeric(Q24), data = dat[dat$pid == "Democrat", ])
summary(lmout.dem.probability)

lmout.io.probability <- lm_robust(probability ~ Trump + Biden + IO + Q21 + Q25 + Q23 + as.numeric(Q24), data = dat[dat$pid == "Independent" |
                                                                                                           dat$pid == "Other", ])
summary(lmout.io.probability)

lmout.rep.probability <- lm_robust(probability ~ Trump + Biden + IO + Q21 + Q25 + Q23 + as.numeric(Q24), data = dat[dat$pid == "Republican", ])
summary(lmout.rep.probability)

terms.probability <- names(coef(lmout.probability))

plot_models(lmout.probability, lmout.io.probability, lmout.dem.probability, lmout.rep.probability,
                    rm.terms = terms.probability[-c(2:4)],
                    show.values = T,
                    show.p = T,
                    legend.title = "Models",
                    m.labels = c("All", "Independents/Others", "Democrats", "Republicans"),
                    vline.color = "lightgrey",
                    p.adjust = "holm",
                    colors = "quadro",
                    digits = 3) +
  labs(title = "") + 
  ylab("Margin of victory: Treatment effects") + 
  xlab("Treatment Effects") + 
  theme_tufte(base_family = "GillSans") 
ggsave("~/Dropbox/curtice-crabtree/plots/ols-margin.png",
       height = 7,
       width = 7)

modelsummary(list(lmout.probability, lmout.io.probability, lmout.dem.probability, lmout.rep.probability),
             output = "latex",
             coef_map = cm,
             stars = T)
